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ABSTRACT 

We present the first vertically resolved hydrodynamic simulations of a laterally propagating, 
deflagrating flame in the thin helium ocean of a rotating accreting neutron star We use a new 
hydrodynamics solver tailored to deal with the large discrepancy in horizontal and vertical 
length scales typical of neutron star oceans, and which filters out sound waves that would 
otherwise limit our timesteps. We find that the flame moves horizontally with velocities of 
order 10^ cm s^^, crossing the ocean in few seconds, broadly consistent with the rise times 
of Type I X-ray bursts. We address the open question of what drives flame propagation, and 
find that heat is transported from burning to unburnt fuel by a combination of top-to-bottom 
conduction and mixing driven by a baroclinic instability. The speed of the flame propagation is 
therefore a sensitive function of the ocean conductivity and spin: we explore this dependence 
for an astrophysically relevant range of parameters and find that in general flame propagation 
is faster for slower rotation and higher conductivity. 
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1 INTRODUCTION 

Type I bursts are tremendous thermonuclear explosions on the sur- 
face of accreting neutron stars (NSs), with luminosities that can 
easily reach the Eddington limit. They are characterized by a fast 
increase in the X-ray luminosity, known as the rise, that lasts from 
less than a few seconds up to ten seconds, and by an exponential, 
slow decay, the tail, that lasts from tens of seconds to few minutes 
( |Lewin et al.|1993[|GalIoway et al.|2008| l. More than ninety known 
X-ray sources have shown Type I bursts (for an up-to-date list see 
In't Zand web page http://www.sron.nl/~jeanz/bursterlist.html). 
All are low mass X-ray binaries, where the NS accretes matter from 
the outer layers of the companion. 

Whether the accreted fluid spreads freely or is confined to part 
of the NS surface probably depends on whether the accretion takes 
place via a boundary layer (the disc directly 'touching' the NS, 
|Inogamov & Sunyaev|2010| l or via magnetic channelling and, sub- 
sequently, on the strength of the magnetic field itself ( Brown &| 
pildstenT998 1. The majority of bursters do not show persistent pul- 
sations, while those that do feature very weak magnetic fields. It is 
therefore reasonable to assume that the accreted material spreads 
over all of the NS surface, forming a thin highly combustible ocean 
consisting of mostly light elements. 

As new fluid piles up, the deeper layers are compressed until 



* E-mail: ycavecchi@uva.nl 



the temperature and density are high enough to trigger nuclear re- 
actions of H, He or both. Depending on the accretion rate and the 
composition of the fluid, the burning can be stable or unstable (jFi>| 
[jimoto et al.|198l] l: in the latter case Type I bursts occur. It seams 
unlikely that the whole star ignites at the very same moment jShara| 
1982); instead it is more likely that the ocean ignites locally and 
that the resulting flame propagates laterally and engulfs the whole 
(or a substantial fraction) of the NS surface. It is the physics of 
the lateral propagation of the thermonuclear flame that is the fo- 
cus of this study. Our goal is to investigate the open question of 
what controls the propagation of the flame and the ignition of un- 
burnt fuel. Mechanisms that may be involved include conduction, 
turbulent mixing associated with convection or other hydrodynam- 
ical instabilities, or compression (of unburnt fuel by elements that 
are already burning and hence expanding). Answering this question 
is critical to explaining burst time scales, such as the rise times. It 
is also relevant to the development of burst oscillations (fluctua- 
tions in the intensity of the burst lightcurves, see |Watts] ( |20I2[ l for a 
review). The mechanisms suggested to explain such oscillations al- 
ways involve some kind of asymmetry on the burning surface, like 
the presence of a hot-spot ( Strohmayer et al. 1996 1 or surface mode 
patterns excited by the flame (Heyl 2004, Piro & B ildsten||2005[ 
[Berkhout & Levin|2008| l. [Spitko vsky et al. (2002, hereafter SLU) 
suggested that the Coriolis force might confine burning material in 
'thermonuclear hurricanes', pointing out the importance of the ro- 
tation of the star and of hydrodynamics for the flame propagation. 
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Early attempts at multi-dimensional simulations of the flame 
propagation mechanism simulations include Nozakura et al 
1984|l, who followed a zonal approach, and|Fryxell & Woosley 



T982ab, who hydrodynamically simulated the first few millisec- 



onds of a detonation in a thick layer. [Zingale et alll ( |2001[ l continued 
this line of research, following the detonating flame for up to hun- 
dreds of milliseconds. However the thick layers required for deto- 
nation are unlikely to accumulate between bursts, so that in reality 
the flame probably develops via deflagration (^Malone et al.|20lT| l. 
The most recent studies of flame propagation, which last again only 
a few milliseconds (far shorter than the timescales of real X-ray 
bursts), are those by |Simonenko et al ( 2012a hi. In the first of these 
papers, the flame propagates via a detonation wave, in a manner 
similar to that in |Zingale et al.|p001f . The authors find that propa- 
gation is due to the hot fluid expanding and spilling over the top of 
the cold fluid, which then compresses and ignites. In the second pa- 
per, the authors explore regimes with densities ~ 2 x lO^g/cm'^ at 
the base of the ocean, conditions where one might expect deflagra- 
tion rather than detonation. However, in this case the simulations 
did not develop a steadily propagating flame. 

None of these studies, however, took into account the rotation 
of the neutron star. SLU studied the role of the Coriolis force on 
flame dynamics and concluded that the Rossby adjustment radius 
may determine the horizontal scale of the burning front. They also 
pointed out the importance of the ageostrophic flow at the hot-cold 
fluid interface as part of the mechanism for the flame propagation 
itself. That said, [SLU] used a two layer, shallow water, method to 
simulate the propagation of the flame and therefore had to make 
phenomenological assumptions about heat and momentum trans- 
port in the vertical direction. Our simulations do not involve such 
assumptions, and are resolved in the vertical direction, allowing us 
to make a detailed study of the flame-propagation physics, taking 
full account of rotation. 

In this paper, we simulate flame propagation on a domain that 
has a horizontal extent that is a substantial fraction of the surface 
of the NS. At the heart of our simulations is the hydrodynamical 
code described in [Braithwaite & Cavecchi|p012[ hereafter BC) . 
It is a multidimensional code, which by construction enforces hy- 
drostatic equilibrium in the vertical direction (this assumption is 
justified since the timescale for sound propagation in the vertical 
direction in the burning layer is much shorter then the nuclear reac- 
tion timescale). Hydrostatic equilibrium allows us to use a longer 
time step, which would otherwise be limited by sound wave prop- 
agation in the vertical direction. It also allows us to use pressure 
as the vertical coordinate. This is a great advantage, because we 
can follow the inflation of the fluid without the need for extra grid 
points which would lie unused for most of the simulation, consum- 
ing extra memory and increasing calculation time (for more details 
see[BCt. 

In Sec. [2] we briefly review the numerical code described in 
and introduce the additions and modifications made to the code 
in order to study thermonuclear flame propagation. We then report 
our results on flame spreading in Sec.|3] We focus in particular on 
the mechanisms that drive flame propagation and investigate the 
speed dependence on rotation rate and conductivity. We conclude 
with a brief summary in Sec.|4] 



2 NUMERICAL IMPLEMENTATION 

In this section we describe the modifications made to the code re- 
ported by BC in order to make it suitable for study of flame propa- 



gation during Type I X-ray bursts. First however we briefly review 
the most salient features of the code as outlined in IBCI It is a 3D 
magnetohydrodynamical code, which uses the cr-coordinate system 
(a pressure coordinate system, see below) on a staggered grid: ther- 
modynamical variables like temperature, pressure, density and heat 
sources are evaluated at the centres of the grid cells, while veloc- 
ities are evaluated on the "faces" of the cells. The code is 3D, but 
for this paper we use a 2D version, assuming variables are indepen- 
dent of one of the horizontal dimensions. We also neglect magnetic 
fields, postponing this for future research. 

The assumption of vertical hydrostatic equilibrium, justified 
by the short vertical sound crossing time (much shorter than any 
time scale of interest in burst simulations), allows us to discard 
vertically propagating sound waves, numerical resolution of which 
consumed the lion's share of the cpu time in previous numerical 
experiments. In our simulations we can therefore employ much 
longer timesteps than previous studies. Vertical equilibrium leads 
naturally to the introduction of a vertical pressure coordinate. This 
in turn makes it possible to follow the fluid as it expands, without 
the need for extra grid cells. 

The code evolves the two horizontal components of the fluid 
velocity, the pressure (which acts as a pseudo density) and the tem- 
perature using a three step Runge-Kutta scheme, while the spatial 
derivatives are calculated with sixth-order finite differences. Pres- 
sure is defined as 



aP, + Pt 



(1) 



where P, = Pb ~ Pt and Pb and Pt are the pressure at the 
bottom and at the top of the simulation. Pt is a constant parameter 
in the simulations and a £ [0, 1] becomes the vertical coordinate 
(0 corresponding with the top, 1 with the bottom). P* can be shown 
to become a pseudo density and the continuity equation becomes 



dP, 

1h 



= -Ia=i 



where 



1= [ V,-(P,u)do-', 



(2) 



(3) 



u being the horizontal component velocity and Vcr = ^dx + ydy, 
taken at constant a. 

Defining the Lagrangian derivative as D/Dt = dt + u^dx + 
Uydy + adcr, the momentum and energy equations are 



DT 1 DP 



VP. 

cr 2nxu + Fvisc. (4) 



(5) 



where 



-gz = P, 



P 



(6) 



p is the density, f2 is the rotation vector of the star, parallel to the 
z direction, Fvisc are the viscous forces, cp is the heat capacity 
at constant pressure and Q is the heat per unit mass per unit time 
(see |BC| for more details); the equation of state in |BC| is a perfect 
monoatomic gas. 

We now move on to discuss the modifications to the code that 
were implemented to render it suitable for Type I burst simulations. 
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2.1 Equation of state 

The first relevant change is to the equation of state (EOS) for the 
fluid in our simulations, which has to be able to describe the physics 
of the NS ocean. We take into account the composition of the fluid 
by expressing it in terms of the mass fraction: X is the fractional 
mass of H, Y that of He and Z = 1 — X — Y the fraction of all 
heavier elements. For a fully ionized perfect gas, the perfect gas 
EOS used in|BC|becomes (assuming full ionization): 



with 



pRT 



with 



12 



7 + 17X + 2Y 



g mol 



(7) 



(8) 



However, since conditions in the NS ocean can lead to elec- 
tron degeneracy, which plays an important role in the vertical sup- 
port of the ocean against the gravitational field, we must take this 
into account in our simulations. We still consider the atoms to be 
fully ionized; however whilst the nuclei are assumed to be a per- 
fect gas, the electrons may be (partially) degenerate and (partially) 
relativistic. We also need to include radiation pressure. 

For this purpose we adapted the publicly available routine 
he Ime o ^of Timmes & S westy 1 2000 1. It uses the density p, tem- 
perature, X and Y to derive pressure, energy, the thermodynamic 
potentials and their derivatives with respect to p, T, X and Y. In our 
code structure p is a derived quantity, while pressure is a primary 
quantity (see |BC| sec. 2). To circumvent the problem of passing 
density as an input parameter to the routine, we interface the orig- 
inal helmeos with a zero-finding routine that calls it repeatedly 
with different values of p until convergence in pressure is achieved. 
Subsequent calls use the previous value of p as an initial guess. 
Given the Courant conditions we impose, the information in each 
grid point does not change much in one time step, and convergence 
is achieved within two or three calls. This is done in parallel for 
each grid point. 

The choice for the EOS is important for the evolution equa- 
tion for the temperature T. We still derive it from the first law of 
thermodynamics: 

DE 



Dt 



p^ Dt 



(where E and Q are the energy and heating rate per unit mass), 
but we have two different results depending on the choice of the 
equation of state. If we use the perfect gas equation of state, the 
evolution equation for temperature has the form (compare to |BC| 
eq. 20): 



Cp DT 
p Dt 



1 DP cpT 
plyt ^ IT 
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.DX 

In 



' Dt ) 



(10) 



where cp = — l)/7, 7 is the adiabatic index and R — 

8.3144621 X lO'' erg K"^ moP ^ is the gas constanQ 

When we include electron degeneracy and radiation pressure, 
by contrast, we have 



^ DT ^DP ^DX ^DY 
''^Dt=''Dt+''-Dt+''l^ 



^ Available at http ://cococubed.asu.edu/code_pages/eos.shtml. 
^ Note that in BC we include p in the definition of R. 
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where we make use of the relations 



(12) 

(13) 
(14) 
(15) 
(16) 

(17) 
(18) 
(19) 
(20) 



since the routine returns variables as functions of p, T, X and Y. 

A further complication comes from the fact that helmeos 
actually uses A = 12/(1 + IIX + 2Y) and Z = A{1 + X)/2, 
(instead of X and Y directly) for P, and that the derivatives of E 
and P are evaluated with respect to T, A, Z and p. Therefore, for 
equations ijTS}, {16}, ^19\ and (|20j we also need 

A^ 
' 12 

' " +(1 + X)^ ) (22) 
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The final evolution equation for the temperature is therefore 
(D/Dt = d/dt + u ■ V,, + dd/da): 



dT/dt = Tt.adv + T .thcrmodyii 

+ Q/cp, 

where the contributions to dT/ dt are separated into 



(u ■ V,r + ddT/da) 



Cp Dt 



Tt,p 



B DX C DY 
5^1jt £^~Dt 



(23) 

(24) 
(25) 

(26) 



and Q/cp, so that we can test the relative importance of the con- 
tributions of the different terms (see Sec. |3.4| (. Q is further divided 
into Q = Qn + Qcond + Qcooi + Qhypcr, where Qn is the nu- 
clear burning contribution, Qcond is the conduction contribution 
and Qcooi is the cooling contribution from the top (see next sec- 
tions). We also include an artificial diffusive term Qhypor (with a 
small coefficient, see Sec. 3.4 of|BC| to ensure numerical stability. 
In the case of the perfect gas EOS the evolution equation is very 
similar to equation l|23j. 

Finally, the term DP/Dt is evaluated according to equation 
(19) of|BC] which does not depend on the EOS, but on the choice of 
the (T-coordinate system. DX/ Dt and DY/ Dt have to be treated 
carefully: in case of reactions, or any change in composition, they 
have to be evaluated explicitly (see Sec.|2.3^. 
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2.2 Conduction 



During burning the composition is evolved according to 



Since conduction may play an important role in flame propagation, 
we include a physical conduction term in Q of the form: 



„d = -V 



p \ 3p/tc 



VT I 



(27) 



where Kc is the effective opacity due to both radiative and conduc- 
tive processes. In the u-coordinate system Qcond takes the form: 



1 



+ 



3k.cP 



dT 



da 



,2dT 

r^r- + — v„ 

oa p 



+ 
2dT 



(28) 



where cf> — gz = P, J^l/pda' . The terms that include (j) are 
due to the fact that the transformation matrix between the Cartesian 
coordinate system and the a coordinate system is not everywhere 
orthogonal. These contributions turn out to be small and can be 
neglected. 

We have implemented two possibilities for evaluating Kc : ei- 
ther it has a fixed value set at the beginning of the simulation, 
or it is calculated for each and every grid point taking into ac- 
count the composition and thermodynamical variables at that po- 
sition. For this second option we adapted the publicly available 
routines of Timmes: sig9^ The opacities calculated in this way 
take into account radiation, scattering and the degree of degeneracy 
(see |Timmes||2000l and references therein). Based on the values 
of density, temperature and composition in the simulations that we 
wanted to perform, we decided to use an average constant value of 
Kc ~ 0.07 cm^ g^^ in our reference simulation, which speeds up 
the calculations whilst still preserving the critical physics. 

As anticipated in section \2A\ we include in equation \23\ a 
hyperdiffusive term (see |BC| section 3.4.2), which mimics conduc- 
tion. This term is unphysical and only used to ensure numerical 
stability. In the horizontal direction, in particular, it will be unphys- 
ically high, and may partly limit the conclusions we can draw from 
our simulations. However test runs involving much lower hyperdif- 
fusivity yielded flame velocities (see section [3^ which differ by 
only a few percent from the values reported in this paper. 
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9y _ 

dt ' 
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dV 



(30) 



(31) 



where the first two terms come from advection and the third is the 
consumption of He due to nuclear reactions {ta = 5.84 x 10^^ 
erg g~^ is the energy production per gram per nucleon). We also 
include a form of artificial diffusion as described in Sec. 3.4 of lBCI 
to ensure numerical stability. 

In terms of sinks of entropy, we include the possibility of cool- 
ing from the uppermost layers. We use a simple formula, derived 
under the assumption that energy is only transported through the 
layers above the simulated computational domain, without addi- 
tional sinks or sources within the atmosphere. This is a somewhat 
coarse approximation, particularly since expansion of the upper 
layers may occur. We use the temperature of the top grid cell in 
order to evaluate the flux due to radiation and conduction: 



p ^ 16(7B ^3 dT 

3pKc dz 



(32) 



where F is the flux, which we assume to be constant in our plane 
parallel approximation and gb and Kc are as defined in Sec. |2.2| 
We further assume that Kc is constant in the layers above the sim- 
ulation. Rearranging and integrating in the vertical, z, direction, 
from the top of the simulation (T) to the top of the NS atmosphere 
(atm), we have 



F pdz 

I atm 



16o-s 



AT . 



(33) 



In hydrostatic equilibrium, the integral on the left hand side reduces 
to Ft / g, so that 



4a-_B y4jT 



(34) 



Then, assuming that the temperature at the top of the atmosphere 
is negligible with respect to that at the top of the simulation, we 
obtain 



2.3 Sources and sinks of heat: nuclear burning and cooling 

Since we are interested in simulating Type I bursts, we implement 
helium burning via the triple-a reaction according to (see , Clayton, 
[1984}: 



= 5.3 X lO^^ps 



-4.4/Tn -1 -1 

e ' erg g s , 



(29) 



where Tg is the temperature in units of 10^ K, Y is the mass frac- 
tion of He and ps is the density in units of 10* g cm"''. Including 
only the triple-a process is of course a simplification, since there 
are many other reaction chains that should be taken into account 
(this model would not be correct even for a pure He accretor). How- 
ever we leave this refinement for later investigation. 



Available at http://cococubed.asu.edu/code_pages/kap.shtml. 
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This is the flux from the surface of a grid cell at the top of the 
simulation. In order to derive the entropy loss per unit mass, we 
multiply the flux by the surface area S of the cell and divide by the 
mass within it: 



Qcoo 



F- 



S 



so that (Az ~ // ~ Pt /gPT) 



Qcool 



PtSAzt 



4g g--B y4 



(36) 



(37) 



This is the sink term we use in our simulations; it could also be 
used as a first approximation to calculate the bolometric lightcurve 
of the bursts. 
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2.4 Tracer particles 

Finally, we add the capability to follow tracer particles. These are 
assigned initial positions uniformly distributed in the integration 
domain and are evolved according to 

^ = Wx(a;,t/,cr) (38) 

^=Uy{x,y,a) (39) 

^=&ix,y,a) (40) 

where d/dt is the rate of change of the particle's position in a coor- 
dinates. Time evolution is the same as for all of the other variables 
(see |BC[ section 3.3), and the values of the three components of 
the velocity at arbitrary points within each grid cell are derived by 
means of bilinear interpolation ( [Press et al.|1992^ of the fluid ve- 
locitjQ 



Run u (Hz) Kc (gem ^) Vf (cms 



I 


450 


1 X 


10+0 


1 -L.vJiJ _i_ U.WOI 




Iq5 


2 


450 


7 X 


10-1 


(1.43 ±0.02) 


X 


10^ 


3 


450 


5 X 


10-1 


(1.52 ± 0.02) 


X 


10^ 


4 


450 


3 X 


10-1 


(1.67 ±0.03) 


X 


10^ 


5 


450 


1 X 


10-1 


(1.91 ± 0.04) 


X 


10^ 


6 


450 


7 X 


10-2 


(2.01 ± 0.05) 


X 


10^ 


7 


450 


5 X 


10-2 


(2.03 ±0.05) 


X 


10^ 


8 


450 


1 X 


10-2 


(1.99 ±0.02) 


X 


10^ 


9 


450 


1 X 


10-^ 


(1.98 ±0.03) 


X 


10^ 


10 


50 


7 X 


10-2 


(1.11 ± 0.18) 


X 


10® 


11 


112.5 


7 X 


10-2 


(5.30 ±0.31) 


X 


10^ 


12 


225 


7 X 


10-2 


(3.04 ±0.10) 


X 


10^ 


13 


900 


7 X 


10-2 


(1.39 ±0.02) 


X 


10^ 



Table 1. Values of the spin frequency u = Q/2n and the opacity Kc used 
in the different simulations. In the third column we report the velocity of 
the flame as measured from the simulations, with en'ors derived from the 
least squares fit (see Section [T2) . 



3 FLAME PROPAGATION SIMULATIONS 

In this Section we describe the numerical setup used for all the sim- 
ulations, then provide a description of what we see in the runs. Fi- 
nally we describe our interpretation of what drives the flame prop- 
agation. 



3.1 Numerical setup 

We ran a series of simulations resolving both the horizontal x and 
vertical z directions, assuming that the dynamical variables are in- 
dependent of the y coordinate (making the simulations effectively 
2D). The fixed initial conditions, common to all of our simulations, 
are: 



Pt 
X 



IQ^^ erg cm 


0.03 



and 



To = 10" K 



Y = 1 
U2 = 0.5 



ST = 3.81 X 10** K 



1) X 10^^ erg cm~^ 



where Pt and P* are the pressure at the top and the difference be- 
tween the bottom and top pressure (see |BC| section 2). Note that 
whilst Pt is constant, P* is a function of horizontal position and 
time, but not of a. The choice of P, means that we simulate 1.7 
scale heights, vi and V2 are the kinetic diffusive coefficients (see 
section 3.4 of |BC[ >. The corresponding coefficients for the temper- 
ature and the composition fractions X and Y are taken to be 1 % of 
these values. 

We also use a common initial temperature perturbation in all 
simulations. We use a ^-independent temperature profile, which 
varies in the horizontal direction according to: 

ST 



To + 



exp[{x - 1.2 kiTi)/0.36 km] 



(41) 



This function ensures that the temperature profile is smooth enough 
that it does not cause numerical issues; 1.2 km corresponds to the 
position where the temperature perturbation of the background To 

■1 We also tested higher order interpolation methods, but found no signifi- 
cant differences. 



is half of its maximum, while 0.36 km is approximately half the 
width between where the perturbation is asymptotic to its maxi- 
mum and where it is asymptotic to its minimum (0 K). 

We simulate a domain with a horizontal extent of 7.5 km, 
which allows more than sufficient room for the propagating flame 
to reach a steady state. The initial conditions have a high tempera- 
ture at one end of the domain, so the flame ignites there and prop- 
agates towards the other end. In some sense the end where ignition 
occurs can be thought of as the 'eye' of the 'anticyclone' system. 
We use symmetric boundary conditions in the vertical direction and 
reflective conditions in the horizontal direction. In all simulations 
presented here, we use horizontal and vertical resolutions of 240 
and 90. Gravitational acceleration g = 2 x W^'^ cm s"2 and we 
use the the plane-parallel approximation and a constant Coriolis 
parameter, i.e. the /-plane approximation. The fluid is at rest at the 
beginning, Ux = cm s-i, and quickly adjusts to the Rossby so- 
lution (see |BC| sec. 4.2) before the flame spreads. 

Since we want to study the effects of different rotation fre- 
quencies and the influence of conduction, we run a series of models 
employing different values of the spin frequency Q and the opac- 
ity Kc. The parameters for the simulations that we run are given in 
Table □ 

To help us find out how the flame propagates, we use test par- 
ticles and follow what happens to these fluid elements before, dur- 
ing, and after ignition conditions are met. We place the test par- 
ticles homogeneously in our grid (note that they are not homoge- 
neous in space since we use a pressure coordinate system) such that 
Xij = i * Sx/200 and cr,^ = j * 1/200, i,j G [1, 200] . We also 
track what happens at three different points in the atmosphere with 
fixed horizontal position. These points rise and descend with time, 
having fixed values of a not z, so this approach is not strictly speak- 
ing Eulerian. However it still allows us to see what happens when 
the flame reaches a determined distance from the ignition point. 



3.2 General description of the propagating flame 

In this section, we give a qualitative description of the burning fluid 
as a whole. We begin by using one particular run as an example, 
since the qualitative behaviour is general. The left hand column 
of Figure [T] shows the fluid in its initial conditions for reference 
run (6). The right hand column shows the conditions at t = 1.15 
s, when the flame is propagating steadily. The upper panels show 
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Figure 1. Initial conditions (left) and conditions at t = 1.15 s (right), when the flame is steadily propagating, for reference simulation (6). Top panels show 
the temperature, middle ones show burning rate with the tracer particles superimposed and the bottom panels show density with isobars superimposed (ten 
levels from P = 10^^ erg cm"^ to P = 6.2 X 10^^ erg cm ''). Note the different horizontal and vertical scales. 
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Figure 2. Electron chemical potential rj at the beginning and at t = 1.15 s for reference run (6). Degeneracy decreases with lower r}. The electrons are always 
partially degenerate, but degeneracy is lifted when the flame passes through. The black line again indicates the position of the inteii'ace. 



the temperature distribution, and the middle ones the burning rate. 
On these panels we superimpose our tracer particles. The bottom 
panels show the density, with isobars superimposed. 

In the top left panel of Figure [T] the fluid is hotter on the left 
of the image: this is the initial perturbation, where the temperature 
is T = 4.81 X 10** K, while at the other side the temperature is 
T = 10* K. Moving to the right panel, we see that the flame front 
has moved to the right. Where the fluid has already burnt, it is hotter 
(T ~ lO" K) and has expanded by a factor of of order fouQ Look- 
ing at the middle panels we can see that the tracer particles have 
been scattered by the passage of the front, while the lower pan- 
els show a drop in density. Eventually, after the flame has passed 
(not shown in the figure), the burning diminishes, the temperature 
decreases and the fluid contracts. 

Looking more closely at the propagating front, we see that 
it is characterized by a slanted interface between the hot burning 
fluid on the left, and the cold unbumt fluid on the right (in Figure 
[Tjright, the interface lies roughly between x = 3.7 x 10^ cm and 
X — 4.7 X 10^ cm). We see a decrease in pressure on the left of 
the interface and an increase immediately to the right of it (see the 
lowest line contours of pressure in the bottom right panel of Figure 
They reflect a change in P* (equation^. Because of the hydro- 
static approximation, is a measure of the column density at each 
point. A change in P* means horizontal mass motion. The decrease 
before the front and the increase after it therefore show that there 
has been a motion of matter from behind the front forward. 

The electrons are partially degenerate, as can be seen from 
Figure [2] where we plot the electron chemical potential. The elec- 
trons remain partially degenerate throughout all the simulations, 
but the degeneracy is lifted by the flame (as can be seen also by 
the fact that the temperature increases by a factor ~ 15, while the 
height of layer increases by only a factor of 4 or 5) so that in the 
hot fluid the perfect gas pressure and the radiation pressure become 
more important. 



^ In general, the maximum expansion factor can be up to ~ 4 — 5 depend- 
ing on the effective opacity Kc that sets the cooling rate. 



The peak of the burning is concentrated in a thin stripe along 
the interface (Figure[T] middle right panel) where the density is still 
high (undiminished by the increase in temperature) and the fuel is 
still almost pure Helium. We also observe tracer particles moving 
in the vertical plane, primarily along the interface and in the region 
to the left of it. 

In order to understand what is driving the flame forward, we 
measure the different terms in the energy equation \23\ : conduc- 
tion Qcond/cp, advection Tt,adv (motion of the fluid) and ther- 
modynamic compression Tt.thermodyn. For each term, we plot the 
contributions in Figure |3] In the four panels the black line approx- 
imates the interface between hot and cold fluid. It is drawn below 
the region of significant burning in order to clearly demarcate re- 
gions where burning has started from those where burrung is about 
to start. 

It is clear from Figure [3] that in the region immediately below 
the peak of the burning, heat conduction is much more important 
in increasing the temperature in the unbumt fuel region than both 
the effects of mixing (measured by the advection of temperature) 
or thermodynamic compression. It is this process that drives flame 
propagation, since the main burning occurs in this zone. In the up- 
per part of the interface, advection and thermodynamic compres- 
sion dominate heat transfer to the unburnt region. That picture is 
confirmed by observing what happens at a fixed horizontal posi- 
tion. In Figure |4] we plot the burning rate versus temperature and 
density, and temperature and density versus time for three different 
positions: at the top, in the middle and at the bottom of the fluid 
at a fixed horizontal position x — 3.3 x 10^ cm. It can be seen 
that the topmost point (black) in the figure is compressed and its 
temperature rises. The lower points then follow, but the burning 
does not really start until the temperature has risen sufficiently. The 
same figure also demonstrates how the burning rate increases with 
increasing temperature, while the correlation with density (see for 
example the lower panel) is not as strong. The decrease of burning 
rate at the end of the curves is due to the consumption of fuel which 
eventually becomes the most important factor. 

Directly above the flame, on the other hand, heat conduction 



8 Y. Cavecchi, A. L. Watts, J. Braithwaite, Y. Levin 















6.15E+11 


9.57E+17 
(k/s) 


1.91E+18 


-9.00E+17 


O.OOE+00 
(k/s) 


9.00E+17 



-9.00E+17 



Qn, t= 1.15E+00S 




(k/s) 

Qcond, t = 1.15E+00S 



















E+17 




O.OOE+00 


9.00 




Qadv, t= 1.15E+00S 




X 10 cm 



-9.00E+17 



O.OOE+00 
(k/s) 



9.00E+17 



Qtherdy, t = 1.15E+00s 



X 10 cm 




X 10 cm 



Figure 3. Burning rate and heating rate associated with advection Tt ^dv conduction Qcond/cp and thermodynamical compression Tt thcrmodyn for 
reference simulation (6). The black line indicates the hot-cold fluid interface. The colour scale has been restricted to highlight details (white regions indicate 
values above the maximum of the scale, black ones values below the minimum). 



is not effective, while the advective and thermodynamic compres- 
sive terms show opposite signs. This is a clear signature of convec- 
tion, which is expected above the burning regions. We note that the 
convective cells near the topmost part of the hot-cold interface are 
mostly parallel to it (i.e. almost horizontal, recall the extreme as- 
pect ratio of the interface), while the ones behind the interface are 
verticaj^ In Figurejsjwe plot contours of the total entropy per unit 
mass as returned by the code helmeos. 

To compute flame propagation speed from our simulations, we 
define the front position as the location with the maximum burning 
rate. In Figure |6] we follow the position of the front for simulation 
(6) and plot it versus time. At the beginning there is a transitional 
stage after the flame is started by the initial perturbation of the tem- 
perature and the front adjusts to a steadily spreading configuration 
(in < 0.1 s). This steady propagation is well fit by a straight line, 

^ We want to stress that also these vertical cells are actually elongated in 
the horizontal direction due to the aspect ratio of our underlying grid cells. 



and the gradient gives us the speed of the flame front. We repeat the 
fit procedure for all of the various runs: the resulting front speeds 
vt are reported in Table^ 

Having measured front velocities, we can determine the ef- 
fects of the rotational spin and the effective opacity Kc (a proxy 
for the heat conductivity). Overall, the gradient of the inflated fluid 
is steeper for higher fi, and the baroclinicity along the interface 
tends to increase with Q.. Flame propagation speed decreases as ro- 
tation rate increases (see next sections and Figure |9]l. Changing Kc 
also has an effect on flame velocity: the flame is faster for lower 
ftc, but the velocity saturates when Kc < 10~^ cm^ g~^ (see Fig- 
ure [TO]l. 

3.3 A first set of conclusions 

Although the fluid moves ageostrophically from behind the inter- 
face forward, this motion does not go past the interface (Figure [T] 
bottom right). We interpret this as follows: the fluid has expanded 
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Figure 4. Top: burning rate versus temperature and density for tliree different points fixed relative to the grid, i.e. at fixed x = 3.3 X 10^ cm and sigma: near 
the top (black), in the middle (blue) and near the bottom (orange). Bottom: temperature and density versus time for the same points. The lines in the lower 
panels have colours corresponding to the same scale as the contours for the burning rate in the middle panel of Figure^ The circles indicate the origin of the 
curves. The strong relation between burning rate and temperature is clear, while the importance of the change of density appears to be less. 



on the left of the front because of its higher temperature, and the re- 
sulting horizontal pressure gradient pushes the hot burning fluid to 
spill over the unburnt fluid. The Coriolis force, however, diverts the 
horizontal x motion into the horizontal y direction and thus creates 
a geostrophic current that compensates for the horizontal pressure 
gradient. The resulting configuration is that of the Rossby adjust- 
ment problem (see |BCp , as anticipated by |SLU| In this case the 
inclination angle a of the interface should be a ~ H/Rb.o, where 
H is the scale height of the fluid and Ruo is the Rossby radius 
^Ro = where Q, = 27rz/ and v is the spin frequency 

of the neutron star. Measuring the slope of the black line in Figure 
[3] we find that the slope is a ~ 3.5 x 10"'', so that its horizontal 
extent is ~ 2 — 3-Rro. This is in accordance with expectations. 

Regarding the motion that we observe in the vertical plane 
along the interface, we note that here the fluid is much more baro- 
clinic, that is to say, the iso-surfaces of density and pressure are 
much more misaligned than elsewhere, as can be seen in the lower 



right panel of Figure[T]and in Figure |7] It is well known from geo- 
physical studies that geostrophic balance is unstable in the presence 
of baroclinicity. The resulting instability is similar in nature to con- 
vection, but with motion, which is no longer vertical, lying within 
the "wedge of instability" between the iso-pressure and iso-density 
contours (|Pedlosky 1 1 987| l. |Fujimoto] ( [T988l \T993\ and |Cumming&l 
|Bildsten| ( |2000 1 in fact already studied the possibility of baroclinic 
instability in the context of Type I bursts, but their baroclinicity 
was very mild since they considered the effects of shear induced by 
the differential rotation due to the vertical expansion of the burning 
layer, and not the effects of the huge horizontal temperature gradi- 
ents that develop during flame propagation. 

In our case, the source of baroclinicity is the horizontally- 
inhomogeneous nuclear burninaj which affects the temperature 



^ Compare Figure^to the middle right panel of Figure^ 
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Figure 5. Entropy per unit mass (radiation, ion and electron gas) for refer- 
ence simulation (6) at t = 1.15 s. The black lines indicate the contours for 
better visualization. The red line indicates the position of the interface. 
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Figure 7. Baroclinicity: VP X Vp, for reference simulation (6). The vector 
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Figure 6. Flame front position for run (6). The symbols indicate the eiTor 
bars on the positions, while the line is the best linear least squares fit. Af- 
ter an initial stage, the flame adjusts to steady propagation. Eventually, the 
flame reaches the opposite boundary (in this case in ~ 3 s). 



profile. Its steady propagation is maintained by the Coriolis force, 
which reinforces the near-geostrophic configuration on time scales 
of order 1/v. Following the tracer particle motion, we can see ad- 
vection along and in front of the interface, which we attribute to 
the development of baroclinic instability. In the previous section 
we noted the presence of cells that are highly elongated in the hor- 
izontal direction at the upper end of the hot-cold fluid interface: 
we identify these cells with baroclinicity-induced motion. FigurejS] 
shows how particles are driven into the front and down along the in- 
terface. After the flame has passed and the front is farther away, the 
particles are caught by the advective motion and driven upwards. 
The ascending part is different for different particles and this pic- 
ture just indicates the general trend. 



ZvsX 

3.4n 



3.1- 




Figure 8. Example of motion of one tracer particle on the vertical plane. The 
red dashed lines indicate the position of the flame front at different times A, 
B, C and D. The con'esponding positions of the particle are indicated by 
the same letters on the particle trajectory. The empty circle indicates the 
starting point. The colours of the trajectory correspond to the same scale as 
the contours for the burning rate in the middle panel of Figure[T] 



3.4 Front propagation mechanism 

Summarizing the results from Section [X2] we see that at the top 
of the interface the fluid is heated up by the spilling over of the 
hot fluid, via advection and thermodynamics. However in the most 
relevant regions for flame propagation, heat is brought across the 
interface primarily by conduction (mainly vertically given the small 
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inclination angl^, with a contribution from baroclinic instability 
mixing. 

The contribution to dT/dt in equation (|23j from conduction 
can be written by means of equation \21) as 



Tt , cond 



t cond 
Cp 



16 gbT 

3 pK, 



3 



vr 



(42) 



From this we can derive an approximate diffusion coefficient for 
conduction -Dcond as 



(43) 



and then derive the timescale for heat diffusion across the vertical 
scale height H as Tcond ~ / Dcond, or 



'7"cond 



Equation (|44j gives 



'^cond 



2.1 X 10"" s 



109 K 



3 p^H^cp 
16 (tbTS ' 



0.07 cm2 g-i 



H 



(44) 



3 X 102 cm 



10^ g cm~^ 
Cp 



108 erg K- 



(45) 



Once the lower fluid has been heated up and starts burning 
it expands again. The Coriolis force then reinforces Rossby ad- 
justment in a time scale of order u^^ ^ Tcond- This translates 
a small vertical shift into a long horizontal displacement, where 
the proportionality is given by the inclination of the interface: 1/a 
~ (2 — 3Rko)/H, as we will see in the next section. 

The effective advective conduction brought about by baro- 
clinic mixing would act on a time scale given by 



Tba 



with 



-D_Lbar ~ 'fXbarAxb 



(46) 



(47) 



where wxbar and Axbar are the physical velocity of the fluid and 
its length scale perpendicular to the hot-cold fluid interface. 

As long as Tbar ^ Tcond, couductiou will be the most effec- 
tive mechanism. In order to get a handle on the importance of the 
baroclinicity induced advection, we measured the average values of 
w±har and Axbar- For each grid point along the interface, we calcu- 
lated the component of fluid velocity in the direction perpendicular 
to the interface. Considering only the region over which wxbar was 
negative, we calculated the average wxbar and measured Axbar as 
the length of this region in the direction perpendicular to the inter- 
face. We then calculated the total average Dxbar: the results are 
reported in Table |2] for all of the simulations. The results confirm 
that baroclinicity is negligible most of the time, apart from in cases 
of very low heat diffusivity (high Kc). One should be aware, how- 
ever, that all these values are order of magnitude estimates, and 
hence only describe trends, not precise timescales. 



8 Heat conduction in the horizontal direction can be neglected since the 
horizontal length scale is larger by a factor ~ 10^ than the vertical length 
scale. Two runs where in one case full conduction was implemented and in 
the other only vertical conduction was used gave virtually identical results. 



Run 


-^Xbar 


'''bar 


'^cond 


''"bar /'^cond 


1 


223.0 


403.6 


142.5 


2.8 


2 


288.6 


311.8 


99.8 


3.1 


3 


328.9 


273.6 


71.3 


3.8 


4 


449.5 


200.2 


42.8 


4.7 


5 


915.0 


984 


14.3 


6.9 


6 


926.0 


97.2 


10.0 


9.7 


7 


1045.1 


86.1 


7.1 


12.1 


8 


2047.2 


44.0 


14 


31.4 


9 


7757.8 


11.6 


0.1 


116.0 


10° 






10.0 




11° 






10.0 




12 


548.4 


164.1 


10.0 


16.4 


13 


1570.9 


57.3 


10.0 


5.7 



Table 2. Diffusion coefficient for the baroclinicity driven advection (equa- 
tion|47| as measured directly from the simulations, its diffusion timescale 
according to equation l|46j and the timescale for conduction as from equa- 
tion j45j, using T = 6 X 10* K and p = 9 X 10^ g cm^^. The last column 
reports the ratio between the two timescales. Note that these values are only 
indicative order of magnitude estimates. 

" These runs had less cleai' configurations, so that reliable measurements 
were not possible. 



3.5 Front propagation speed 



Following e.g. [Landau & Lifshitz| ( |1959| >; [Fryxell & Woosieyl 
l |1982b[ l, the velocity of flame propagation across the interface 
should be given, in a deflagration regime, by 



VfA 



(48) 



where Tn is the burning time scale, given by — ta/Qn, see 
Equations \29\ and ^0[, 



1.1 X 10"^y"^exp(4.4 X lO'' K/T) 



VlOSKy V gem 



(49) 



In order to estimate the horizontal propagation velocity across the 
NS, this velocity has to be multiplied by the factor (2 — 3)i?Ro/-ff 
which expresses the ratio of the area of the burning front to the area 
of the vertical section of the ocean (see [Landau & Lifshitz|1959l l. 
The horizontal velocity becomes 
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9 1 


3cp p'^Tn 





vi ~ 1.8 X loV^''^ exp(-2.2 x lO'"* K/T) cm s"^ 
V450 Hz/ 



0.07 cm2 g-i 



H 



3 X 102 cm 



-1/2 



2 X 1014 cm s-2 

1/2 



Cp 



108 ergg-1 K-i 



(50) 



1/2 



(51) 



Again, if Tbar ~ Tcond, thcu the actual v[± should be given by 
a combination of conduction and advection, with an extra term of 
order 



(52) 



to be included in equation l [48[ ). The result should then be multiplied 

by the same factor, (2 — 3)Rko/H. 
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Figure 9. Velocity of flame propagation versus l/v. All these runs have 
Kc = 0.07 cm2 g-i. 
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Figure 10. Velocity of flame propagation versus l/y^K^. All these runs 
have V = 450 Hz. 

The scaling expected from equation \5l) agrees with what we 
measure in our simulations. The agreement is within half an order 
of magnitude and this allows us to put constraints on the numerical 
factors in front of equation ijSTJ which are not determined by the 
order of magnitude estimates that led to it. 

We verify the dependence of v{ on spin rate u using runs 6-10. 
In Figure|9] we can see that increasing the rotation frequency slows 
down the flame, as expected, with a.l/v dependence. 

By contrast Figure [TO] where we plot Vi against the inverse 
heat conduction (our effective opacity Kc). shows more complex 
behaviour. For opacities Kc 0.05 g cm~^, the flame speed in- 
creases approximately with 1 / as expected from equation \5i) . 
Below it, the velocity seems to asymptote to Vf ~ 1.99 x 10^ cm 
s~^. In our simulations we see a change in the morphology of the 
flame: indeed all simulations showed a flame leaning on the hot- 
cold fluid interface similar to the middle right panel of Figure [T] 
apart from simulation (9) (Kc ~ 0.001 g cm~^), which does not 
show such a leaning flame and has a much more vertical structure of 
the temperature profile (Figure [TT]l, and simulation (8) (kc = 0.01 
g cm^^) where this trend is beginning to become apparent. 

We interpret this point as marking the transition where the 
conduction time scale Tcond becomes comparable to the burning 
time scale Tn and the thickness of the slanted burning front becomes 
comparable to the vertical scale height. At this point the front speed 
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Figure 11. Temperature and reaction rate for simulation (9) after the flame 
is steadily propagating. The morphology of the flame is different to that 
shown in Figure[T] This simulation shows the asymptotic behaviour seen in 
Figure [To] 



saturates and the whole layer burns through on the timescale Tn: 
according to the original estimate ofjSLUj(see the discussion lead- 
ing to equations 21-23), this means that the horizontal speed sat- 
urates at ~ -Rro/tii. If we adopt the values T = 6 x 10* K and 
p = 9 X 10^ g cm^^, as in Table[2] we obtain an average burn- 
ing rate of Tn ~ 0.45 s, so that 2 Rro/tu ~ 1.9 x 10^ cm s~\ 
This order-of-magnitude estimate is in good agreement with what 
we measure. 

On the other hand, when Kc ^ 0.3 g cm~^, we observe a 
greater deviation from our estimates. In this limit of small conduc- 
tivity, we suspect that the baroclinic motions become more impor- 
tant. Their contributions, on top of those predicted by equation \5l\ , 
have to be taken into account, until the front speed asymptotes to a 
baroclinical-motion driven system. 
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4 DISCUSSION AND CONCLUSIONS 

In this paper, we have been able to simulate for the first time the lat- 
eral propagation of a deflagrating vertically resolved flame on the 
surface of a neutron star. We find that after an initial post-ignition 
adjustment, the front propagates steadily with constant velocity, un- 
til it reaches the opposite side of the simulation box. The fact that 
the flame velocity is constant (Figure |6| gives us confidence that, 
regardless of the physics of localised ignition, steady flame prop- 
agation depends only on the physics acting in the ocean layer and 
the conditions therein. After all the surface has been traversed by 
the flame, the fluid columns cools down slowly, in a time which 
depends on the opacity (see Equation |37[(, whilst still burning the 
residual fuel. We note that in 2003, Anatoly Spitkovsky (unpub- 
lished) obtained somewhat similar flame fronts using the Pencil 
code. Due to computational constraints, however, he assumed un- 
physically large neutron star spins, so that the Rossby radius was 
comparable to the ocean scale height. The micro-physics of the 
flame propagation mechanism was not identified, and full explo- 
ration of parameter range was not carried out (Spitkovsky, private 
communication) . 

We have explored the dependence of the flame speed on the 
spin frequency of the neutron star v and the heat conductivity of 
the fluid (expressed as an inverse of the effective opacity Kc). We 
measured velocities in the range 1.33 x 10^ - 1.11 x lO'' cm s~^, 
which cross the entire domain of 7.5 km in 0.7 - 5.6 s. These num- 
bers are in good agreement with the rise times observed from Type 
I burst sources, suggesting that we have included all the relevant 
physics in our simulations and that we are now in a position to 
explore in more detail the behaviour of flame propagation during 
Type I bursts. 

The flame propagates through a combination of the 
ageostrophic forward flow of the burning fluid on top of the as-yet 
unbumt fluid (as argued previously in SLU I, and top-to-bottom heat 
transport across the large-area strongly-inclined interface between 
burning and cold fluid. Heat transport leading to ignition is effected 
primarily by microscopic heat conduction and, in runs where con- 
ductivity was set to lower values, by baroclinic motions. 

In Section [33] we derived an order of magnitude estimate for 
the velocity that the front would have if it were driven by conduc- 
tion. We calculated a dependence of the speed on Kc and v of the 
form 1 / v^/k^ (Figures|9]and 10 1 and confirmed these expectations 
with the results of our simulations. A breakdown of this Kc depen- 
dence is seen at both low and high Kc, which can be understood 
qualitatively. In particular, we observe the existence of a possible 
asymptote in the velocity when the effective opacity is too small, 
which we explain as follows. When the opacity decreases suffi- 
ciently, the conduction time scale becomes shorter than the nuclear 
burning time scale. The latter becomes the bottleneck, the burn- 
ing front width becomes comparable to the scale height, and the 
nuclear burning time scale becomes the time scale of vertical ex- 
pansion. This translates into a horizontal velocity of ~ i?Ro/T"n, as 
already anticipated by |SLU| 

There are a number of hydrodynamical issues that now have 
to be explored further. Firstly, the effect of the baroclinic instability 
at the hot-cold fluid interface could be explored in more detail. 
Secondly, the flow in the j/-direction has a velocity comparable to 
the sound speed, and Kelvin-Helmholtz instabilities that might be 
generated by this flow need to be investigated. Other aspects of the 
flame propagation will be explored in future work, including the 
effects of a better burning prescription taking into account elements 
other than Helium. We also aim to investigate the possibility of 



exciting large-scale waves in the ocean, and the effect of magnetic 
fields. 

Finally, some of our simulations suggest that in the absence 
of a sufficiently strong Coriolis force the flame will die out. This 
leads to an important question: can the flame cross the equator? 
Near the equatorial belt the effective Coriolis force is much weaker 
and this could lead to rapid lateral spreading of the burning front, 
and subsequent quenching of the burning by enhanced cooling. 
This would have important consequences for efforts to determine 
NS radius from observations of type-I bursts (see, e.g., |Steiner| 
|et al.|[2010^ , since it is usually assumed that the whole star is 
burning at the peak, and the derived radius of the burning area is 
used as a measure of the NS radius. If the flame cannot cross the 
equator, this fact has to be taken into account when dealing with 
those estimates. This would also have important implications for 
burst recurrence times, and may help to explain the properties of 
multi-peak bursts ( [Bhattacharyya & Strohm ayer| |2006"| . We plan 
to investigate this possibility by introducing a variable Coriolis 
parameter in future work, to simulate properly the changes that 
would occur as a flame approaches the equatorial belt. 
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